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Lunar globules are smooth glassy objects which were discovered by the astronauts on the 
moon. These objects are small - most are less than a tenth of a millimeter in diameter - 
though some are a good deal larger. The shapes of the globules vary on a continuum from 
spheres to prolate spheroids to dumbbells. Figure 1 contains two photographs, one of a 
prolate globule and one of a dumbbell shaped globule. 

Several theories have been proposed for their creation. All these theories assume that 
molten rock is shot from the surface of the moon, solidifies in space above the moon and 
then falls back to the surface. 

The rotational theory which we study in this paper makes the following assumptions: The 

volume of the molten rock does not change during cooling. The angular momentum is con- 
served. There are no internal motions because of the high viscosity of the molten rock, 
i.e., in equilibrium the globule is rotating as a rigid body. Finally, we assume that the 
kinetic reaction of the globule to the forces is fast relative to the rate of cooling, 
i.e., we assume that the globule reaches equilibrium at constant energy. For more discus- 
sion of this theory see [1]. 

Most thermodynamists would agree that with these assumptions equilibrium will be 
achieved at a local minimum of the internal energy and that the relevant part of the inter- 
nal energy is 


E = jl (n)u 2 + T • AOS1) 


Here n is the region occupied by the globule in a coordinate system fixed in the body. 
I(n) is the moment of inertia of the globule with rest>ect to its axis of rotation. A(3Q) 
is the surface area of the globule. T is a positive constant which describes the surface 
tension properties of the molten rock (this number is effectively constant over wide varia- 
tions of temperature) , and to is the angular velocity. 

We remark that from the standpoint of rational thermokinetics there is a need to justify 
these claims theoretically. However, we will assume here that equilibrium will be achieved 
at a local minimum of this internal energy expression. 

We therefore consider the mathematical problem: 


Minimize : 

E = i Iu 2 + T • A Oft) 

(i) 

Subject to: 

Volume (n) = V(ft) = k^ * given constant 



I u * k 2 * given constant. 


sometimes more 

convenient to consider the problem in the 

following equivalent form. 
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Minimize: 

E = — + AOfi) 

(2) 

Subject to: 

V(n) * k^ = given constant 



The class of globules over which we minimize will always be some subset of the smoothly 
bounded axisyrametric (with respect to the axis of rotation) configurations. Here axisyra- 
metric with respect to the z-axis means that if (x,y,z) is a point in the globule then so 
is (-x,-y , z) . 


We would first like to point out that for this variational problem there is no configu- 
ration for which E (or E) attains its absolute minimum, so that we can only expect to find 
a local minimum. To see this, consider problem (2). Let e > 0 be given. By detaching 
two small identical pieces from the sphere at volume k^ and placing them symmetrically on 
either side of the remainder of the sphere on an axis which is perpendicular to the axis of 
rotation, it is possible by pushing the two small pieces far away from the sphere to obtain 

a configuration, Q , such that E(q) ■ ai + e where ai is the surface area of the sphere of 

volume k^. Thus the absolute minimum for £ is ai if it exists. But clearly no configu- 
ration has this value for E provided that k 2 ¥ 0, since che sphere of volume kj is the only 

possibility. 

Applying the classical methods of the calculus of variations we obtain as the Euler- 
Lagrange equation for this problem the equation: 

C (p) - - (Ar 2 + B) f (3) 

where C(p) is the mean curvature of the surface 3n at the point p e , r 2 * x 2 + y 2 where 
we have taken the z-axis is the axis of rotation, and A and B are constants which involve 
Lagrange multipliers and must be chosen to satisfy the constraints. It is reassuring to 
note that this equation can be obtained directly by balancing the forces at each point on 
3 ft without considering the variational problem at all. Since the pressure inside the 
globules must be greater than the outside pressure, we may conclude that the constants A 
and B are positive. Here we have adopted the convention that inward curvature is negative. 

It is possible to write this equation as a free boundary value problem for a nonlinear 
partial differential equation. But this problem seems difficult to treat even numerically. 
We will introduce additional symmetries into the problem so that we are led to consider an 
approximate problem involving ordinary differential equations. 

The oblate solutions 


There exist exact solutions to (3) which are surfaces of resolution with respect to the 
z-axis (the axis of rotation). By taking advantage of the assumed symmetries, we see that 
these solutions will be solutions to the free boundary value problem, 

f f" = (1 + (f) 2 )(l - 2f (1 + (f) 2 ) 1 * (Af 2 + B)) (4) 

f ' (0) = 0 : f (b) = 0, f ' (b) = - 

Here, f and b are unknown and A and B are positive constants which must be chosen to sat- 
isfy the constraints. The function x * f(z) generates the surface of the globule when 
rotated about the 2-axis. Remarkably, this equation can be integrated exactly in terms of 
elliptic integrals. (See Chandrasehkar [2].) However, it is easier to do numerically. 
These solutions while they are exact solutions to the Euler-Lagrange equation (3) , they 
almost surely are not local minima for our variational problem except possibly near the 
sphere i.e., for small values of the angular momentum. Nevertheless, there is a smooth one 
parameter family of these solutions for a given volume starting with the sphere (k? = 0) 
and becoming more and more oblate. A convenient parameter to index these solutions is, I, 
their moment of inertia. Table 1 contains some important numbers for a few members of this 
family . 


The approximate prolate and dumbbell solutions 

Motivated by the fact that the actually occurring shapes are close to being surfaces of 
revolution with respect to an axis (which we call the x-axis) which is perpendicular to the 
axis of rotation, we consider the variational problem modified so as to include in the 
class of possible globules only those which are surfaces of revolution with respect to the 
x-axis. Now the appropriate functional to minimize is 

rb 

J ^ f 4 * 2 + f?x ' w? + 2n T f (f ) ? 


+ x ^ nf 2 + Xj t^- f k + pit f 2 x 2 ]iii dx 

where the function y * f (x) generates the surface of the globule by rotation about the 
x-axis, p is the density, and X] and \2 are Lagrange Multipliers. Here o> is treated as a 
parameter. As before, both f and b are unknown. This leads to the free boundary value 
problem, 
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f f" 


(5) 


- (1 + (f ’ 7 )(1 - 2f(l + (f’ 7 )‘* (Ax 7 ♦ jf } + B)' 


f • (0) -Os f (b) - 0, f • (b) - - 

As before, the positive constants A and B must be chosen to satisfy the constraints. As in 
the oblate case, if we fix the volume, we can generate, numerically, a continuous one par- 
ameter family of solutions starting at this sphere and proceeding through prolate shapes to 
dumbbell shapes with narrower and narrower necks. Again, I, the moment of inertia is an 
increasing parameter along this family of deformations of the sphere. Figure 2 contains 
graphs of some of the members of this family. The volumes of all the globules in this 
graph are equal. The reader is asked to compare the shapes in this figure with the photo- 
graphs in Figure 1. Mathematically, this equation is more difficult to handle than ( 4 ), 

It cannot be integrated in terms of elliptic integrals and no existence theorem is known 
for this free boundary v*r.lue problem. Table 2 contains extensive information about this 
family. The number e is defined as the thickness at the axis of rotation divided by the 

length of the globule. The number, D, is defined by D * - ~ . We note that both of these 

B 3 

numbers are invariant under similarity transformations and that if a globule is a solution 
to our free boundary value problem then so is any similar globule. Both Tables 1 and 2 
were computed assuming that the volume is one and that the density is one. 

We have, therefore, two families of deformat ions of the sphere each parameterized by, I, 
the moment of inertia of the globule with respect to its axis of rotation . We denote the 
oblate family by Uq{T) and the prolate family by up(I), For given I, the two corresponding 
members are local minimum for the variational problem modified in two ways. First, we add 
the constraint that I « const. Second, we restrict the variations allowed to be such that 
they produce in the case of the oblate family surfaces which are surfaces of revolution 
with respect to the axis of rotation (the x-axis) and in the prolate family surfaces which 
are surfaces of revolution with respect to the x-axis. We caution that if more general 
variations are allowed the prolate family globule cannot be a local minimum since they do 
not satisfy (3). While the oblate family globule is also not a local minimum, at least 
sufficiently far from the sphere, even though it does satisfy (3). 

The first constraint, I * constant, is eliminated by considering E as a function of I 
along either of the families. We have 

ckv 

E (I ) * + A ( T ) . 


At a minimum for the problem where the constraint, I * constant, is dropped but the 
restrictions on the variations are retained we must have -jy - 0. That is ! ? jjy * c 

d * F 

■Mso we must have Some simple calculations show that this condition is equivalent 

to the condition (I'gy) n 0 assuming that yy » 0 at the point considered. An examina- 
tion of Table 2 leads to the conclusion that these globules are no longer a minimum beyond 
qlobule #25. That is beyond this point the globules are clearly unstable. If this theory 
is correct, the dumbbell shaped globule in Figure 1 is near the limit of stability. If it 
were much narrower at the neck it would break apart. 

Compa rise n of t he energies 

We assume now that this volume of all globules is fixed at one. If k2 is sufficiently 
small, there is exactly one local minimum for the modified variational problems along each 

of the families Uq(I) and up (I). The value of E at a point where gy * 0 is given by the 

expression Tjy + A which forms the last column in the tables. By inspecting Tables 1 and 2 

we see that for a given k 2 the energy, Ep(k2) of the minimizing globule from the flp(I) 
family is smaller than the energy Eo(k 2 ) at the minimizing globule from the Agd) family 
from globules #16 or #17 on to the end of the table, while up to that point Eq^) is 
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smaller than Ep(kp), However, since there exist purturbat ions of any of the prolate 
family globules which produce an even smaller energy there are globules near the prolate 
family which produce smaller values of E nearer to the sphere than globule #16. 

We conjecture that there is a one parameter family of globules which start from the 
sphere and develops essentially like the prolate family but which contains true local 
minimum for the variational problem (1). (Chandrasehkar (21 expects that the oblate 
spheroids are stable near the sphere and that a bifurcation occurs along this family which 
contains ellipsoids with different lengths for all three axes.) 

Con e I us ion 

The qood agreement between the actual shapes of the Umar globules and the numerical 
results in Figure 2 certainly lend some support to the rotational theory. However, since 
the existence of true local minima has not been established and since the theoretical 
understanding of the thermokonetics is still to be attained, there is certainly room for 
doubt . 


Th e nu meric al o omp ut at ion 

The simplif icat ion of the variational problem yields a second order ordinary differen- 
tial equation (ODE) for f (x) 


d-f 

,Tx- 


1 

fTxT 


(1 ♦ 


•S' 


') 1 - 2f(x) / 


1 + 


',ar; 

( ax’ 


(Ax- + B * 0 <f (x) ) 




(6) 


and one must find a positive value of the independent variable, say Xf, so that the fol- 
l ow i ng bou rid a r y co nd i t ions ho l d 


d t 


x*0 


0 ; 


f(x f ) - 0; 


df i 
dx 


x=x. 


(?) 


i.e., the function should start off flat from the f-axis and should cross the x-axis going 
straight down. 

The numerical solutions present in this paper were computed using 14 digit precision 
BASIC on a 2-80 based microcomputer employing two routines from the small machine oriented 
library of mathematical software, SCRUNCH [4] to solve this variant of a two-point boundary 
value problem by simple shooting. 7F.R01N, a robust root finder developed by L. F. Shampine 
and R. C, Allen, Jr. which uses a careful combination of bisection and the secant rule, was 
used to find the missing end point, Xf. RKF4S, a fourth/fifth order Runqe-Kutta-Fehlberq 
method originally coded by H. A. Watts and h. F. shampine for solving the initial value 
problem for systems of first order differential equations with automatic step selection 
and reliable 'rror control, integrated the DDEs for each trial value of Xf. 

The differential equations and boundary conditions were reformulated in terms of arc 
lenqth and integrated from x*xf to s*0 using a suggestion of C. W. dear to avoid the sin- 
gularities present at the boundaries and the instability due to a large, positive eigen- 
value of the ODE system (6) at x*xf. Changing the direction of integration caused the sys- 
tem to be initially stiff due to the large, neqativo eigenvalue present, but this causes no 
problem due to the small step sizes typically necessary to start the integration. Lettinq 
x » x(s); y « y(s) • f(x(s)), the problem actually computed was 


d : x 

Ss? 


(Ax' + B + Cy J ) 


dy ,dx 1 , 

at ( 3? yTST 



,dx 

2 cfi 


(Ax' ♦ B + 


cy') 


dx 

3i 


,dx 1 . 

( 3S yTsT 


( 8 ) 


with boundary conditions at x f ■ x(0): 
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x(0) - x f , g 


s-0 


0; y (0) - 0; 


S-0 


( 9 ) 


and defining the total length of the arc, s^, by the condition x(s T ) - 0 


x(s T ) » 0; 


dx 

5? 


x«s n 


-1? y(s T ) - positive? 


s»s n 


0. 


The singularity still present at the starting point x f due to the term a? yhr in (8) ia 
avoided by using the approximation 


( 10 ) 


/dx 1 > 
( 3? yTsT ) 


* - (Ax 2 + B + Cy 2 ) 

s«0 x=x f 

y*o 


(Ax 2 + B) 


which comes from using the first two terms of a Taylor series for ^ and y(s) near s*0 
and the known properties of the problem. 

The final BASIC program takes as input the values of the parameters A and B (using the 
relation that C = A/8) and an interval (xi,X2) in which to search for the unknown Xf, the 
root of the function 




s=s„ 


obtained by integrating (8) backwards from the trial starting point x using initial condi- 
tions (9) until the function crosses the y-axis. The differential equations were inte- 
grated using a mixed relative and absolute error tolerance of 10~ 6 and the root was found 
to a similar mixed tolerance of 10“ 5 . These tolerances were easily met using the 14 digit 
BASIC now that the problem is formulated in a numerically stable manner. Once the missing 
value Xf was obtained, the initial value problem (8 and 9) was solved and the values (x,y) 
obtained at each step of the integration were plotted using a simple line-printer/plotter 
routine. 


Table 1 


c 

V-l 

20 

25 

50 

75 

150 

I 

.12768 

1.9094 

.19978 

.20477 

.21300 

A 

4.8833 

4.9297 

4.9742 

5.0022 

5.0523 

dA 

dl 

3.499 

5.034 

5.611 

6.087 

7.107 


.1105 

.1835 

.2239 

.2553 

,3235 

x ar + A 

5.505 

5.891 

6.065 

6.249 

6.566 
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Table 2 


1 

e 

D 

I 

A 

dA 

ar 

12 ar 

T dA 4. A 

x ai + A 

1 

1 

0 

.15401 

4.8363 

.26 



2 

1 .05 

.125 

.15554 

4.8367 

1.09 



3 

1.07 

.25 

.15691 

4.8382 

1.98 



4 

| 

1.10 

.5 

15934 

4.8430 

2.59 



5 

1.14 

.75 

.16181 

4.8494 

3.22 



6 

1.21 

1.25 

.16694 

4.8659 

4.14 



7 

1.43 

2.5 

.18516 

4.9413 

4.52 



8 

1.45 

2.625 

.18815 

4.9548 

4.55 



9 

1.47 

2.6375 

.18848 

4.9563 

4.57 



10 

1.48 

2.6875 

.18986 

4.9626 

4.46 

.161 

5.81 

11 

1.49 

2.71875 

.19078 

4.9667 

4.54 

.165 

5.83 

12 

1.50 

2.75 

.19175 

4.9711 

4.50 

.166 

5.83 

13 

1.55 

3.05 

.20779 

5.0433 

4,38 

.189 

5.95 

14 

1.92 

3.0 

.23741 

5.1729 

4.13 

.233 

6.15 

15 

2.02 

2,875 

.25021 

5.2257 

3.99 

.250 

6.22 

16 

2.15 

2.6875 

.26628 

5.2898 

3.87 

.275 

6.32 

17 

2.15 

2.6375 

.26770 

5.2953 

3.88 

.278 

6.33 

■>ITT73»T^>e WIM 

2.18 

2.625 

.27139 

5.3096 

3.78 

.278 

6.34 

19 

2.26 

2 . 5 

2.8150 

5.3478 

3.62 

.287 

6.37 

20 

■EH 

2.25 

.30217 

5.4227 

3 42 

.312 

6.46 

21 

2.58 

2.0 

.32432 

5.4984 

3.20 

.336 

6.54 

22 

2.67 

1.75 

3.4908 

5.5775 

2.98 

.363 

6.62 

23 

2.95 

1.5 

.37758 

5.6623 

1 

2.75 

.392 

6.70 

24 

3.21 

1.25 

.41136 

5.7552 

2.34 

.396 

6.72 

25 

3.98 

.75 

.50533 

5,0755 

1.98 

.506 

6.98 

26 

4.66 

.5 

.57291 

6.1059 

1.52 

.500 

6.98 

27 

6.30 

.25 

.65168 

6.2258 

1.02 

.431 

6.88 

28 

9.15 

.125 

.68487 

6.2595 
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Figure 1 

Photographs supplied courtesy of Dr. Gertrude Hinsch, Department of Bioloqy, University of 
South Florida and Dr* H* Fertig, Max Plank Institute fur Kerm Fhysik, Heidelberg, Germany* 



Figure 2 

The numbers appearing on the graph are values of D. 
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